clear; clc; close all;

% 参数同之前
q = 0.75;  
a = 2;     
P = 2;
Q = 1;
x0 = 2;
y0 = -2;

b_min = -5;
b_max = 5;
M_b = 100;
b_list = linspace(b_min, b_max, M_b);

N_total = 1500;
N_transient = 500;

FraAppEn_list = zeros(1, M_b);

% 计算过程
for idx = 1:M_b
    b = b_list(idx);
    [x, y] = SCLMM(q, a, b, P, Q, x0, y0, N_total);
    y_vec = y(N_transient+1:end);
    
    % 计算分数阶近似熵
    FraAppEn_list(idx) = fractional_approx_entropy(y_vec, 2, 0.1*std(y_vec));
    
    fprintf('b=%.3f, FraAppEn=%.5f\n', b, FraAppEn_list(idx));
end

% 绘制图像
figure('Position',[200 200 900 600]);
plot(b_list, FraAppEn_list, 'g-', 'LineWidth', 1.5);
xlabel('b', 'FontSize', 14);
ylabel('Fractional Approximate Entropy', 'FontSize', 14);
title(sprintf('Fig.5(c) 2D-SCLMM 分数阶近似熵 (q=%.2f)', q), 'FontSize', 16);
grid on;
box on;

%% 分数阶SCLMM函数
function [x, y] = SCLMM(q, a, b, P, Q, x0, y0, N)
    x = zeros(1, N);
    y = zeros(1, N);
    x(1) = x0;
    y(1) = y0;

    w = zeros(1, N);
    for k = 1:N
        w(k) = exp(gammaln(k + q - 1) - gammaln(k));
    end

    delta_x = zeros(1, N);
    delta_y = zeros(1, N);

    for n = 2:N
        delta_x(n-1) = a * sin(P * y(n-1)) * sin(P * x(n-1)) - x(n-1);
        delta_y(n-1) = b * cos(Q * x(n-1)) * cos(Q * y(n-1)) - y(n-1);

        sum_x = 0;
        sum_y = 0;
        for i = 1:n-1
            sum_x = sum_x + w(n - i) * delta_x(i);
            sum_y = sum_y + w(n - i) * delta_y(i);
        end

        x(n) = x0 + sum_x / gamma(q);
        y(n) = y0 + sum_y / gamma(q);
    end
end

%% 近似熵函数 (m为嵌入维数, r为容忍阈值)
function ApEn = fractional_approx_entropy(U, m, r)
    N = length(U);
    
    % 构建嵌入维m和m+1的序列块
    Xm = zeros(N-m+1, m);
    Xm1 = zeros(N-m, m+1);
    
    for i = 1:(N-m+1)
        Xm(i,:) = U(i:i+m-1);
    end
    for i = 1:(N-m)
        Xm1(i,:) = U(i:i+m);
    end
    
    % 计算距离并计数
    Cm = zeros(1, N-m+1);
    Cm1 = zeros(1, N-m);
    
    for i = 1:length(Cm)
        dist = max(abs(Xm - Xm(i,:)), [], 2);
        Cm(i) = sum(dist <= r) / (N - m + 1);
    end
    for i = 1:length(Cm1)
        dist1 = max(abs(Xm1 - Xm1(i,:)), [], 2);
        Cm1(i) = sum(dist1 <= r) / (N - m);
    end
    
    % 计算phi
    phi_m = sum(log(Cm)) / (N - m + 1);
    phi_m1 = sum(log(Cm1)) / (N - m);
    
    % 近似熵
    ApEn = phi_m - phi_m1;
end
